Anomalous specific heat of supercooled water 
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Abstract 

We study by Monte Carlo simulations and mean field calculations the specific heat Cp of a 



(N 



-C cell model that reproduces known properties of water. In the supercooled regime we find two Cp 

maxima that are strikingly similar to very recent experimental results. We show that the maximum 
at higher temperature is associated with the fluctuations of hydrogen bond formation and that the 

a 

second maximum at lower temperature is associated with the cooperativity of the hydrogen bonds. 
If this cooperativity is not taken into account, the lower temperature maximum disappears. Our 
finding offers a thermodynamic interpretation for the two maxima in Cp. 
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Because of their in physics, chemistry and biology, the anomalies of water are objects of 



intense research 



BE 



2]. One of water's anomalies is its large isobaric specific heat Cp, that 



a. 



increases upon cooling below 35°C 

Onepossible interpretation of the anomalies of water is the singularity free (SF) sce- 
nario 4|, according to which water's anomalies are a thermodynamic consequence of a 
negatively sloped locus of temperatures of maximum density in the pressure-temperature 
(P-T) plane. In this scenario Cp has a maximum that does not change value with P but 
shifts to lower T. 

A second scenario is the hypothesized existence of a first order transition line separating 
a low density liquid (LDL) and a high density liquid (HDL). By moving along this line, 
the density difference between LDL and HDL decreases and disappears at a liquid-liquid 
critical point (LLCP) (P c , T c ) [5J. The loci of maxima of thermodynamic response functions 
converge as the LLCP is approached to the Widom line (the locus of maxima of correlation 
length) [6j. 



A third scenario 
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called the critical-point-free (CPF) scenario, hypothesizes the 
presence of an "order-disorder" first order phase transition, without a critical point. 

Recent experiments on water confined in cylindrical silica gel pores with diameters of 
1.2-1.8 nanometers allow one to probe extremely low temperatures that are inaccessible 
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to bulk water. One finds two maxima in Cp as the temperature decreases [9|, [l0|, 

prominent peak at low T is accompanied by a smaller and broader peak at higher T. These 

I I 

experiments have been interpreted in terms of non- equilibrium dynamics [U]. Here we offer 
an alternative thermodynamic interpretation of this phenomenon, based on Monte Carlo 
(MC) simulations and analytic mean-field (MF) calculations for a cell model [12] that is 
able to reproduce the three scenarios described above for water [l3J. Our interpretation is 
supported by very recent experiments [14 ] . 

The system consists of N particles in a <i-dimensional volume V, which is divided into 
N equivalent cells each with one molecule i G [1 , N] and with volume V/N larger than a 
hard-core volume vq. The interaction Hamiltonian is 15 
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The first term with J > accounts for the hydrogen bond energy, where Oij = 1, . . .,q 
are the Potts variables representing the four bond indices of molecule i with respect to 



four nearest-neighbor (n.n.) molecules j, 5 a ^ = 1 if a = b and 5 a ,b = otherwise, and 
< i, j > denotes that i and j are n.n. The second term represents an intramolecular (IM) 
interaction accounting for the 0-0-0 correlation [l7[ , locally driving the molecules toward a 
tetrahedral configuration. When the bond indices of a molecule are in a state corresponding 
to a local tetrahedral arrangement, the energy is decreased by an amount J a ^ 0, and (k, l) t 
indicates one of the six different pairs of the four bond indices of molecule %. The third term 
Uw{r) denotes the van der Waals interaction between molecules at distance r = (V/N) 1 ^, 
represented by a Lennard- Jones potential with attractive energy e > J and truncated at a 
hard-core distance r = (vq) 1 ^. 

Experiments show that the formation of a hydrogen bond leads to a local volume expan- 
sion [if]. Thus in our system the total volume is V = Vmc + N hb vhb, where Vmc ^ Nv is 
a dynamical variable allowed to fluctuate in the simulations, and 

Nhb = J2 <W^ ( 2 ) 

is the total number of hydrogen bonds, vhb is the constant specific volume increase due to 
the hydrogen bond formation. 

We perform MC simulations for the two-dimensional case, with N = 10 4 molecules, each 
in a cell with four n.n. cells, at constant P and T. The MC dynamics consists in updating the 



variables a^ by means of the Wolff algorithm 18) and updating the volume V in accordance 



with the acceptance probability min (l,exp [— (3 (AE + PAV — NkBTln(Vf/Vi))]), where 
AE is the variation of the right hand side of Eq.flT]) with the update, (3 = (fc^T) -1 , fee is the 
Boltzmann constant, Vi and Vf are respectively the initial and final values of the volume, 
and AV = V f - Vi. We find a LLCP at P c v /e = 0.70 ± 0.02 and k B T c /e = 0.053 ± 0.001 
for the choice of parameters q = 6, J/e = 0.5 and J a /e = 0.05. We study pressures in the 
interval 0.001 ^ Pv /e ^ 1.5. 

We calculate Cp = (dH/dT) p , where H = (E) + P(V) is the enthalpy, and (•) denotes 
the thermodynamic average. For low pressure isobars, such as Pvo/e = 0.001, we observe 
the presence of two Cp maxima: one, at higher T, and the second, at lower T, sharper 
[Fig. [U(a)]. The less sharp maximum moves to lower T and eventually merges with the 
sharper maximum as P is raised toward P c . The temperature of the sharper maximum does 
not change much with P at low P; its value slowly increases, reaching the largest values 
at the critical pressure P c [19]. Approaching P c from below the two maxima merge. For 



P > P c this maximum occurs at the temperature of the first-order liquid-liquid (LL) phase 
transition. For P ^> P c the two maxima split: Cp for the sharper maximum decreases in 



value and shifts to lower T along the LL phase transition line, whi 



e Cp for the less sharp 



maximum is independent of P [Fig. QJb)], as has been noted 20l. |21|. 

We also calculate Cp in the MF approximation [7|.We find that the two maxima are 
distinct only well below P c [Fig. [2(a)] . Both maxima move to lower T as P increases, though 
the less sharp maximum at higher T has a more pronounced P-dependence. Above Pv^/e ~ 
0.3, the two maxima merge into a single maximum. We also find that for higher P [Fig.[2](b)] 
the maximum of Cp increases on approaching the MF critical pressure P c MF i>o/e = 0.81±0.04 



and that the single maximum for P > p^' IF marks the LL phase transition line 19|, |22 |. 

To understand the origin of the two Cp maxima, we write the enthalpy as the sum of two 
terms 

H = H SF + H Coop (3) 

where H SF = (-JN HB + P{V M c + N HB v HB )) and H Coop = H- H SF . Hence, we consider 
C P = C S P F + Cp° op where we define the SF component C| F = (dH SF /dT) p and the 
cooperative component C p oop = (dH Coop / dT) [Fig. 131(a)]. Cf, F is responsible for the broad 
maximum at higher T. Cp F captures the enthalpy fluctuations due to the hydrogen bond 
formation given by the terms proportional to the hydrogen bond number Nhb- This term 
is present also in the SF model [4J . To show that this maximum is due to the fluctuations of 
hydrogen bond formation, we calculate the locus of maximum fluctuation of Nhb, related to 
the maximum of \dNnB/dT\p [Fig. IU(a)], and find that the temperatures of these maxima 
correlate very well with the locus of maxima of Cp F [Fig. [5] . 

We find in Fig. [3](a) that the maximum of Cp at lower T is given by the maximum 
of C P oop . To show that C P ° op corresponds to the enthalpy fluctuations due to the IM 



term in Eqfj] proportional to J a |23j, we calculate \dN IM /dT\ P , where N IM is the number of 
molecules with complete tetrahedral order. We find that the locus of maxima of \dNiM/dT\p 
[Fig. IU(b)] overlaps with the locus of maxima of C p oop [Fig. |3j. Therefore, the maximum of 
C p oop occurs where the correlation length associated with the tetrahedral order is maximum, 



a. 



i.e. along the Widom line associated with the LL phase transition 

In MF we may compare Cp calculated for the LLCP scenario ( J a > 0) with Cp calculated 
for the SF scenario ( J CT = 0) [Fig. [31(b)]. We see that the sharper maximum is present only in 
the LLCP scenario, while the less sharp maximum occurs at the same T in both scenarios. 



We conclude that the sharper maximum is due to the fluctuations of the tetrahedral order, 
critical at the LLCP, while the less sharp maximum is due to fluctuations in bond formation. 



The similarity of our results with the experiments in nanopores is striking 11] . Data 



in ll|] show two maxima in Cp. They have been interpreted as an out-of-equilibrium dy- 
namic effect in [2|, QjJ, but more recent experiments 14( show that they are a feature of 



equilibrated confined water. Therefore, our interpretation of the two maxima is of consid- 
erable interest. 

We thank F. Mallamace, S. Sastry and E. Strekalova for discussions, NSF grant 
CHE0616489 and Spanish MEC grant FIS2007-61433 for support. 
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FIG. 1: (a) Temperature dependence of the specific heat Cp from MC simulations, for the param- 
eters in the text, along low pressure isobars with P < P c . A broad maximum is visible along with 
a more pronounced one at lower T. The first maximum moves to lower T as the pressure is raised 
and it merges with the low-T maximum at Pvo/e ~ 0.4. Upon approaching P c v$/e = 0.70 ± 0.02 
the sharp maximum increases in value, (b) Same for P > P c : the two maxima are separated only 
for Pvo/e > 0.88; the sharp maximum decreases as P increases. In both panels errors are smaller 
than symbol size. 
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FIG. 2: Same as in Fig. [T] but from mean-field calculations (a) at P < P* and (b) at P > P ( 
The mean-field critical pressure is P^ IF = 0.81 ± 0.04. 




FIG. 3: (a) Decomposition of Cp from MC simulations [Fig. [I] for Pvq/e = 0.1 into the cooperative 
component C p oop and the SF component Cp F . (b) Comparison of MF calculations for the LLCP 
scenario case (J CT /e = 0.05) and the SF case (J CT = 0). The low-T maximum is present only in the 
LLCP case. Both lines are calculated at Pvo/e = 0.1 
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FIG. 4: (a) Temperature dependence of (|dA^^s/dT|)p for different isobars, (b) Temperature 
dependence of (\dNiN/dT\)p for different isobars. 



10 



1.5 r^ 



c 
> 

Oh 




FIG. 5: Phase diagram from MC simulations showing the liquid-gas transition, the liquid-liquid 
transition and the temperature of maximum density (TMD). Emanating from the LLCP is the 
Widom line, the locus of maxima of Cp F , locus of maxima of dNns/dT and dNjM/dT. 
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